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ABSTRACT 

Electron  tomography  allows  determination  of  the  three- 
dimensional  structures  of  cells  and  tissues  at  resolutions  sig¬ 
nificantly  higher  than  is  possible  with  optical  microscopy. 
Electron  tomograms  contain,  in  principle,  vast  amounts  of 
information  on  the  locations  and  architectures  of  large  num¬ 
bers  of  subcellular  assemblies  and  organelles.  The  develop¬ 
ment  of  reliable  quantitative  approaches  for  interpretation 
of  features  in  tomograms,  is  an  important  problem,  but  is 
a  challenging  prospect  because  of  the  low  signal-to-noise 
ratios  that  are  inherent  to  biological  electron  microscopic 
images.  As  a  first  step  in  this  direction,  we  report  meth¬ 
ods  for  the  automated  statistical  analysis  of  HIV  particles 
and  selected  cellular  compartments  in  electron  tomograms 
recorded  from  fixed,  plastic-embedded  sections  derived  from 
HIV-infected  human  macrophages.  Individual  features  in 
the  tomogram  are  segmented  using  a  novel,  robust  algo¬ 
rithm  that  finds  their  boundaries  as  global  minimal  surfaces 
in  a  metric  space  defined  by  image  features.  Our  expecta¬ 
tion  is  that  such  methods  will  provide  tools  for  semi- automated 
detection  and  statistical  evaluation  of  HIV  particles  at  differ¬ 
ent  stages  of  assembly  in  the  cells,  and  present  opportunities 
for  correlation  with  biochemical  markers  of  HIV  infection. 

1.  INTRODUCTION 

Transmission  electron  microscopes  have  conventionally  been 
used  in  biomedical  research  to  obtain  two-dimensional  pro¬ 
jection  images  of  thin  objects  such  as  molecules,  cells  and 
tissues.  Such  images  can  be  recorded  in  most  modern  elec¬ 
tron  microscopes  at  magnifications  ranging  from  ~100x  to 
~l,000,000x.  The  use  of  electron  microscopes,  is,  how¬ 
ever,  not  limited  to  imaging  in  2D.  Using  emerging  methods 
in  electron  tomography  (see  [1]  for  a  recent  review),  it  is 
now  also  possible  to  routinely  determine  three-dimensional 
structures  using  principles  that  are  very  similar  to  those  used 
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in  technologies  such  as  computerized  axial  tomography.  Thus, 
one  can  record  a  series  of  images  of  a  given  object  over  a 
wide  range  of  tilt  angles,  and  combine  them  using  back  pro¬ 
jection  algorithms  to  generate  three-dimensional  volume  of 
the  imaged  object. 

A  key  problem  in  biological  electron  tomography  is  that 
the  images  obtained  are  at  relatively  low  signal-to-noise  ra¬ 
tios.  In  part,  this  is  because  of  the  tremendous  complexity 
of  biological  specimens;  for  example  a  single  human  cell 
can  contain  thousands  of  copies  of  tens  of  thousands  of  pro¬ 
teins  packaged  in  a  variety  of  multi-protein  complexes  and 
organelles  of  differing  shapes  and  sizes.  A  second  factor 
comes  from  the  potential  of  electrons  to  damage  organic 
matter,  which  necessitates  the  use  of  electron  doses  that  are 
high  enough  to  obtain  measurable  contrast,  but  low  enough 
to  minimize  structural  damage.  The  rapid,  quantitative  in¬ 
terpretation  of  the  vast  amount  of  data  in  tomograms  of  cells 
and  tissues  therefore  poses  a  challenging  problem.  This 
issue  is  becoming  increasingly  important  now  because  of 
the  rapid  advances  in  instrument  automation  that  have  led 
to  dramatic  enhancements  in  the  speed  of  data  collection. 
We  are  interested  in  developing  approaches  for  3D  segmen¬ 
tation  of  features  in  cellular  tomograms  that  can  work  ro¬ 
bustly  and  rapidly  despite  the  low  signal-to-noise  ratios.  As 
a  test  case,  we  have  used  tomograms  recorded  from  human 
macrophages  infected  with  HIV.  The  cells  were  fixed,  em¬ 
bedded  in  plastic,  stained  with  uranyl  acetate  and  lead  cit¬ 
rate  and  sectioned  in  an  ultramicrotome  to  produce  sections 
with  thicknesses  in  the  range  of  150nm  to  200  nm.  These 
sections  were  placed  on  an  electron  microscopic  grid  coated 
with  a  thin  carbon  film,  and  imaged  in  a  Tecnai  12  electron 
microscope  operating  at  120  kV  equipped  with  a  LaB6  fil¬ 
ament.  Tomograms  were  constructed  using  standard  back- 
projection  algorithms  as  implemented  in  the  IMOD  recon¬ 
struction  package  [2]. 

Figure  1  shows  slices  from  a  tomogram  recorded  from 
a  small  region  of  cells  infected  with  HIV.  Within  this  slice, 
there  are  several  identifiable  features  which  bear  a  resem¬ 
blance  to  the  slice  of  either  an  assembled  virion  or  enclosed 
membranous  entities  with  varying  interior  density  relative 
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to  the  cytoplasmic  medium.  Our  goal  is  to  detect  these 
structures  with  minimal  user  bias,  analyze  them,  and  estab¬ 
lish  correlation  of  the  nature  and  extent  of  these  features 
with  progression  of  viral  infection.  In  this  work,  we  con¬ 
centrate  primarily  on  the  segmentation  of  these  features  and 
basic  statistical  analysis  of  their  distribution  in  the  volume. 


Fig.  1.  Four  types  of  features  identified  manually  in  a  slice 
extracted  from  a  tomogram.  The  tomogram  was  obtained 
from  a  chemically  fixed,  plastic-embedded,  stained  200  nm 
thick  section  from  an  HIV-infected  macrophage.  Each  of 
them  has  an  envelope  that  is  stained  more  darkly  than  the 
background.  Some  such  as  the  one  pointed  to  by  the  solid 
arrow  resemble  an  assembled  HIV  particle,  while  others, 
such  as  those  pointed  to  by  the  dashed  arrows  are  vesicular 
in  nature  and  could  represent  entities  that  are  either  precur¬ 
sors  or  packaging  units  for  the  fully  assembled  virus. 


2.  3D  TOMOGRAPH  SEGMENTATION 


to  eliminate  irrelevant  objects  that  are  not  meaningful  for 
the  subsequent  analysis.  We  should  also  note  that  at  the 
present  stage  of  the  investigation,  our  goal  is  to  select  as 
many  features  as  possible  in  order  to  minimize  user  bias  in 
particle  selection. 

2.1.  Segmentation  of  individual  features 

The  segmentation  problem  is  formulated  as  the  search  for 
a  minimal  surface  on  a  metric  space  that  depends  on  image 
features,  with  the  restriction  that  it  contains  the  specified 
interior  point.  Let’s  S  denote  the  3D  surface  representing 
the  boundary.  The  minimal  surface  problem  is  stated  as  [3]: 

S  =  argmin  /  g(Sp)  ds  (1) 

Jsp 

where  the  metric  g  :  5ft3  [0,  oo]  will  signal  whether  or 
not  a  point  is  in  the  boundary  of  a  feature  element.  Specif¬ 
ically,  g  will  take  low  values  for  points  located  over  thin 
dark  regions  and  higher  values  elsewhere  in  the  volume. 
The  minimal  surface  S  is  then  encouraged  to  go  through 
areas  of  small  g ,  avoiding  regions  of  high  g.  As  in  every  en¬ 
ergy  minimization  setting,  two  issues  have  to  be  addressed: 
first,  find  the  appropriate  g  metric  so  the  surface  will  follow 
the  features  of  interest;  and  second,  solve  the  minimization 
problem  to  get  the  globally  optimal  surface. 


As  pointed  out  above,  electron  tomographic  images  gener¬ 
ally  display  very  low  signal-to-noise  ratios,  and  the  task  of 
segmenting  features  of  interest  can  be  a  challenging  one. 
We  are  developing  new  algorithms  that  can  work  under  these 
conditions,  as  described  below. 

As  all  the  features  of  interest  have  a  dark  boundary  around 
them  (i.e.  only  the  inside  is  different),  we  will  use  the  same 
segmentation  algorithm  to  get  the  boundary,  and  at  the  sec¬ 
ond  stage,  classify  them  by  looking  at  the  interior  voxels. 
Based  on  very  simple  criteria  (i.e.  average  gray  value)  we 
can  separate  filled  from  empty  structures,  and  more  com¬ 
plex  criteria  can  be  used  to  achieve  a  finer  classification  (see 
Section  4). 

The  tomograms  are  segmented  in  a  semi-automatic  fash¬ 
ion:  a  single  point  inside  each  cell  is  first  specified  by  the 
user  selecting  only  those  structures  that  correspond  to  sim¬ 
ple  closed  boundaries  (this  is  done  by  inspection  of  all  slices 
in  the  3D  volume).  For  each  selected  point  we  segment 
the  surrounding  structure  using  a  robust  segmentation  al¬ 
gorithm  described  next.  Although  the  detection  of  interior 
points  can  be  done  automatically  (e.g.,  from  singularities 
of  the  distance  function  to  robust  edges),  the  structure  of 
the  tomograms  is  so  complex  (membranes  merged  together, 
boundaries  with  large  gaps,  presence  of  secondary  structural 
elements,  etc.),  that  it  will  require  intensive  post-processing 


2.1.1.  Designing  the  image  dependent  metric 

The  goal  is  to  define  a  metric  that  captures  the  local  geom¬ 
etry  of  the  image  and  signals  the  presence  of  flat  dark  areas 
which  correspond  to  3D  boundary  voxels.  Following  the 
work  in  [4]  we  look  at  eigenvectors  and  eigenvalues  of  the 
Hessian  matrix  to  characterize  the  local  structure  of  an  im¬ 
age.  Let’s  Xk  denote  the  eigenvalue  with  the  k- th  smallest 
magnitude  (i.e.  |Ai|  <  | A2 1  <  | A3 1).  Plate-like  structures 
have  one  predominant  eigenvalue  and  corresponding  eigen¬ 
vector  in  the  direction  normal  to  the  plane  (the  direction 
where  gray  value  changes  the  most).  The  other  two  eigen¬ 
values  have  vanishing  values  and  the  corresponding  eigen¬ 
vectors  form  an  orthogonal  basis  of  the  plane.  By  looking  at 
the  sign  of  the  first  eigenvalue,  we  can  tell  apart  dark  plate¬ 
like  features  (Ai  >  0)  from  bright  ones  (Ai  <  0).  Com¬ 
bining  these  factors  together  we  define  a  geometric  measure 
function  A4  :  5ft3  -A  [0,  00]  as: 

U  =  !  0  Ai  >  0 

I  v/ATTAfTA|-|Ai|  otherwise 

In  a  post-processing  step,  a  non-maxima  suppression  algo¬ 
rithm  keeps  only  local  maxima  of  M  in  the  direction  per¬ 
pendicular  to  the  plane  (this  is  the  3D  extension  of  the  tech¬ 
nique  used  in  the  Canny  edge  detector,  see  [5]).  By  keeping 


only  local  maxima  in  the  direction  normal  to  the  plane  we 
obtain  a  better  localized  feature  indication  function. 

Observe  that  the  bigger  the  value  of  M  at  a  given  point, 
the  stronger  the  indication  that  the  point  may  be  on  the 
boundary  of  a  structure.  We  then  choose  the  metric  g  to 
be  g(M)  =  (1  -  Mn)  +  w,  where  Mn  =  ma%M)  is  the 
normalized  measure  (0  <  Mn  <  1),  and  w  is  a  small  con¬ 
stant  that  prevents  g  from  vanishing,  see  [6].  As  required, 
g  1  +  w  for  background  voxels  (Mn  0)  and  g  — ^  w 
for  points  located  on  the  boundaries  of  interest  (Mn  — 1). 

2.7.2.  Finding  the  minimal  surface 

Once  the  metric  is  defined,  we  are  ready  to  solve  the  mini¬ 
mization  problem  and  get  the  surface  S  that  minimizes  the 
energy  (1).  This  problem  is  usually  solved  in  a  variational 
framework,  we  compute  the  Euler-Lagrange  equations  and 
find  the  optimal  way  to  reduce  the  energy  (1)  given  an  ini¬ 
tial  energy  state.  Unfortunately,  this  will  only  find  a  local 
minimizer  (the  closest  one  to  the  initial  condition)  and  can¬ 
not  guarantee  the  convergence  to  the  global  minima.  Many 
algorithms  have  been  proposed  in  the  literature  that  pro¬ 
vide  different  mechanisms  trying  to  drive  the  surface  to¬ 
wards  the  global  minima,  e.g.  [7,  8].  Although  they  offer 
improved  performance,  they  can  still  get  trapped  in  local 
minima.  There  are  a  few  approaches  that  guarantee  conver¬ 
gence  to  the  global  minima  in  the  2D  case  (curves  on  2D 
images),  see  [6,  9].  For  the  3D  case,  the  recent  work  [10] 
finds  the  global  minima  for  surfaces  but  requires  a  pair  of 
curves  that  lay  on  the  object  as  initialization. 

Our  approach  for  finding  the  minimal  surface  is  inspired 
by  [9],  extended  to  surfaces  instead  of  planar  curves.  The 
only  input  required  by  the  algorithm  is  a  point  inside  the 
cell.  The  surface  minimization  problem  is  then  solved  im¬ 
posing  the  restriction  that  the  initial  point  should  be  inside 
the  surface.  A  3D  polar  coordinate  transformation  is  first 
performed  with  center  in  the  point  inside  the  cell.  We  then 
look  for  an  initial  surface  guess  in  polar  space,  using  a  tech¬ 
nique  similar  to  the  one  in  [11].  The  resulting  surface  (that 
we  express  in  implicit  representation)  is  back-converted  from 
polar  coordinates  to  the  Cartesian  grid  along  with  the  in¬ 
trinsic  distances  (computed  in  the  polar  domain  with  the  g- 
metric).  Back  on  the  Cartesian  grid,  we  denote  f  the  im¬ 
plicit  surface  representation  and  d  the  back  converted  dis¬ 
tances.  As  a  refinement  step,  f  is  evolved  according  to 
the  following  partial  differential  equation:  <pt  =  k,  \V(j)\  + 
Mn  (Vd  •  V0),  where  Mn  is  the  normalized  feature  indi¬ 
cator  function  defined  before.  The  curvature  term  is  a  regu¬ 
larity  constraint  and  the  external  advection  field  Vd  pushes 
the  surface  towards  the  minima  of  the  energy  (1).  The  fac¬ 
tor  Mn  turns  the  advection  term  off  in  voxels  that  are  not  on 
the  boundary  (Mn  0),  so  pure  regularization  will  nicely 
fill-in  the  missing  gaps. 


The  segmentation  for  each  3D  structure  is  run  indepen¬ 
dently  and  in  a  serial  fashion.  The  approximate  running 
time  for  each  of  them  is  less  than  3  minutes  in  a  1.2  Mhz 
laptop  computer.  The  processing  of  the  whole  volume  can 
well  be  parallelized  to  significantly  reduce  the  computation 
time.  Some  examples  are  shown  in  the  next  section. 

3.  RESULTS 

In  Figure  2  we  show  representative  resulting  3D  surfaces  su¬ 
perposed  on  top  of  slices  of  an  unprocessed  tomogram.  The 
segmented  surface  shows  an  excellent  fit  to  the  boundaries 
of  the  vesicular  features.  We  also  show  that  the  segmenta¬ 
tion  algorithm  can  be  used  to  classify  the  volumes  in  terms 
of  the  mean  internal  density,  as  illustrated  in  a  segmented 
2D  slice,  where  the  regions  are  automatically  classified  (red 
and  green)  based  on  differing  internal  average  grey  values. 
We  also  show  a  2D  slice  of  the  full  tomogram  with  the  seg¬ 
mentation  curves  superposed,  Figure  3. 


Fig.  2.  Two  examples  of  3D  segmentation  of  cell  structure 
generated  from  a  point  inside  specified  by  the  user.  The 
right  example  shows  volumes  corresponding  to  low  density 
(green)  and  to  high  density  (red). 


Fig.  3.  We  show  a  slice  of  the  3D  volume  with  the  segmen¬ 
tation  curves  on  top.  Two  density  classes  are  shown  (red 
and  green)  automatically  classified  according  to  the  average 
gray  level  inside  the  volumes. 

Once  all  relevant  structures  are  segmented,  we  can  per¬ 
form  some  simple  statistical  analysis  on  the  results.  As  each 
volume  is  obtained  as  an  implicit  function,  geometry  com¬ 
putations  (e.g.  size,  average  gray  value,  shape,  etc.)  are  eas¬ 
ily  obtained.  In  Figure  4  we  show  histograms  of  the  average 


gray  level  (density)  distributions  inside  the  selected  volumes 
in  two  different  tomograms.  Note  that  we  can  clearly  clas¬ 
sify  cells  into  two  classes:  the  ones  with  the  filled  interior 
and  the  empty  ones.  Furthermore,  looking  at  the  spatial  dis¬ 
tribution  of  gray  values  inside  each  cell,  more  sophisticated 
criteria  can  be  devised  to  classify  in  the  different  cell  types 
presented  in  Section  1 . 


Fig.  4.  Top:  Distribution  of  average  gray  levels  inside  seg¬ 
mented  regions  in  two  different  tomograms  from  different 
regions  of  the  infected  cell.  Bottom:  Distribution  of  size 
in  segmented  volumes  within  a  tomogram,  normalized  with 
reference  to  the  largest  volume  in  the  set. 

Figure  4  also  shows  the  size  distribution  of  cells  within 
each  tomogram.  Note  that  as  overall  thickness  of  the  tomo¬ 
gram  is  only  ~200  nm,  and  the  virion/vesicular  entities  are 
~100  nm  wide,  only  a  few  cells  are  captured  completely  in¬ 
side  the  volume.  Acquisition  of  serial  tomograms  from  suc¬ 
cessive  sections  that  are  stitched  together  computationally 
should  allow  analysis  of  larger  volumes,  therefore  provid¬ 
ing  more  reliable  statistics  on  the  segmented  volumes. 

4.  CONCLUSION 

In  this  paper  we  have  report  a  new  algorithm  specifically 
designed  for  the  semi- automated  segmentation  of  features 
in  cellular  tomograms  obtained  by  electron  microscopy.  We 
demonstrate  this  is  applicable  for  the  identification  and  pos¬ 
sible  classification  of  HIV  particles  in  infected  macrophages, 
and  could  provide  a  quantitative  basis  for  analyzing  HIV  in¬ 
fection. 

We  are  currently  acquiring  new  data  to  study  the  evo¬ 
lution  of  the  statistics  (see  for  example  Figure  4)  when  the 
virus  evolves.  In  addition  to  simple  statistics  such  as  aver¬ 
age  gray  value  (representing  the  density  inside  the  region) 


and  volume,  we  plan  to  carry  out  shape  analysis  and  classi¬ 
fication  of  the  segmented  regions  using  novel  computational 
approaches  for  shape  statistics  being  developed  in  the  com¬ 
puter  vision  literature. 
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